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A multilayered particle is illuminated by plane acoustic or electromag- 
netic waves of one or several frequencies. We consider the inverse scat- 
tering problem for the identification of the layers and of the refraction 
coefficients of the scatterer in a non-Born region of scattering. Local de- 
terministic and global probabilistic minimization methods are studied. A 
special Reduction Procedure is introduced to reduce the dimensionality of 
the minimization space. Deep's and the Multilevel Single-Linkage methods 
for global minimization are used for the solution of the inverse problem. 
Their performance is analyzed for various multilayer configurations. 
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1. INTRODUCTION 

Many practical problems require an identification of the internal structure of an 
object given some measurements on its surface. In this paper we study such an 
identification for a multilayered particle illuminated by acoustic or electromagnetic 
plane waves. Thus the problem discussed here is an inverse scattering problem. A 
similar problem for the particle identification from the light scattering data is stud- 
ied in [29]. The precise formulation of the problem is postponed till Section 2. Our 
approach is to reduce the inverse problem to the best fit to data multidimensional 
minimization. This is done in Section 3. It is also shown there that more than one 
frequency of the incoming waves is required to provide a stable identification. The 
resulting minimization is a challenging problem, since the objective function has 
many narrow local minima. Finding a global minimum (the sought identification) 
is the main subject of the study here. In Section 4 we analyze various local mini- 
mization methods and develop a special Local Minimization Method. This method, 
together with a specially designed Reduction Procedure, is capable of finding this 
type of local minima. In Section 5 Rinnooy Kan and Timmer's Multilevel Single- 
Linkage Method for global minimization is presented. It is paired with the Local 
Minimization Method of Section 4, and, finally, gives the tool for the successful 
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scatterer's identification. A detailed numerical evidence of the performance of this 
method is presented in Section 6. 

2. DIRECT PROBLEM 

Let D C R 2 be the circle of a radius R > 0, 

D m = {x e R 2 : r m _! < \x\ < r m , m = 1, 2, . . . , N} (2.1) 

and S m — {x e K 2 : = r m } for = r < ri < • • • < r^ < R. Suppose 
that a multilayered scatterer in D has a constant refractive index n m in the region 
D m , m = 1, 2, . . . , N. If the scatterer is illuminated by a plane harmonic wave 
then, after the time dependency is eliminated, the total field u(x) = Ui{x) + u s (x) 
satisfies the Helmholtz equation 

Au + fcp«i = 0, \x\ > r N (2.2) 

where Ui(x) — e tk " x ' a is the incident field and a is the unit vector in the direction of 
propagation. The scattered field u s is required to satisfy the Sommerfeld radiation 
condition at infinity, see [8]. 

Let k m = fc 2 n m . We consider the following transmission problem 

Au m + k m u m = x eD m , (2.3) 

under the assumption that the fields u m and their normal derivatives are contin- 
uous across the boundaries S m , m = 1,2, . . . ,N. 

In fact, the choice of the boundary conditions on the boundaries S m depends on 
the physical model under the consideration. The above model may or may not be 
adequate for an electromagnetic or acoustic scattering, since the model may require 
additional parameters (such as the mass density and the compressibility) to be ac- 
counted for. However, since the goal of this paper is to study algorithms capable to 
resolve the Inverse Scattering Problem, we will accept the above simplified problem 
here. For more details on transmission problems, including the questions on the 
existence and the uniqueness of the solutions, see [27], [1] and [13]. 

The Inverse Problem to be solved is: 

IPS: Given u(x) for all x G S = {x : \x\ = R) at a fixed fc > 0, find the 
number N of the layers, the location of the layers, and their refractive indices 
n m , m = 1,2, ... ,N in (2.3). 

Here IPS stands for a Single frequency Inverse Problem. Numerical experience 
shows that there are some practical difficulties in the successful resolution of the 
IPS even when no noise is present. While there are some results on the uniqueness 
for the IPS (see [1]), assuming that the refractive indices are known, and only the 
layers are to be identified, no stability estimates are available. The identification 
is successful, however, if the scatterer is subjected to a probe with plane waves of 
several frequencies. Thus we state the Multifrequency Inverse Problem: 

IPM: Given u p (x) for all x e S — {x : \x\ = R) at a finite number P of wave 
numbers > 0, find the number N of the layers, the location of the layers, and 
their refractive indices n m , m = 1, 2, . . . , N in (2.3). 
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3. BEST FIT PROFILES 

If the refractive indices sufficiently close to 1 , then we say that the scat- 

tering is weak. In this case the scattering is adequately described by the Born 
approximation, and there are methods for the solution of the above Inverse Prob- 
lems. See [8], [9], [23] and [24] for further details. However, when such an assumption 
is inappropriate, the preferable method is to match the given observations to a set 
of solutions for the Direct Problem. Since our interest is in the solution of the IPS 
and IPM in the non-Born region of scattering, we choose to follow the best fit to 
data approach. This approach is used widely in a variety of applied problems, see 
e. g. [4]. 

Note, that, by the assumption, the scatterer has the rotational symmetry. Thus 
we only need to know the data for one direction of the incident plane wave. For 
this reason we fix a = in (2.2) and assume that the (complex) data functions 



gW(0), p=l,2,...,P (3.1) 

are given for < 9 < 2ir, corresponding to the observations measured on the 
surface S of the ball D for a finite set of free space wave numbers kjf^ . 
Fix a positive integer M. Given a configuration 



Q = {ri,r 2 , ■ ■ ■,r M ,n 1 ,n 2 , ■ ■ . ,n M ) (3.2) 

we solve the Direct Problem (2.2)-(2.3) (for each free space wave number k^) 
with the layers D m = {x e R 2 : r m _i < |x| < r m , m — 1, 2, . . . , M}, and the 
corresponding refractive indices n m , where r = 0. Let 

w^(e)=u^(x)\ xes . (3.3) 

Fix a set of angles 9 = (81,62, . ■ . , 0l) and let 



HI2 = (5> 2 (*o) 1/2 ( 3 - 4 ) 

1=1 

Define 

1 P || — g( p ) || 2 

$(ri,r 2 , . . • ,r M ,n 1 ,n 2 , . . . ,n M ) = -p \\g(p)\\ 2 — ~' 

where the same set Q is used for as for w^ p ' . 

We solve the IPM by minimizing the above best fit to data functional $ over an 
appropriate set of admissible parameters A a d m C K 2M . 

It is reasonable to assume that the underlying physical problem gives some esti- 
mate for the bounds ni ow and rihigh of the refractive indices n m as well as for the 
bound M of the expected number of layers N. Thus, 
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Aadm C {(n,r 2 , ■ ■ ■ ,r M ,n 1 ,n 2 , ■ ■ ■ ,n M ) ■ < n < R , ni ow < n m < n high } . 

(3.6) 

Note, that the admissible configurations must also satisfy 

ri < r 2 < r 3 < ■ ■ ■ < r M ■ (3.7) 

As it was already mentioned in Section 2, the numerical evidence shows that IPS 
is, practically, unresolvable. Here is an example to illustrate the situation. Let 
the configuration Q 1 be (0.4,0.6,0.49,9.0) with N = 2 and R = 1.0. Thus Q 1 
corresponds to the two layer cylinder 



n(x) 



0.49 < x < 0.4 
9.0 0.4 < |^| < 0.6 
1.0 0.6 < Id < 1.0 



Let Q 2 = (0.3794,0.5662,0.6377,0.040,8.282,5.969) with N 
thus Q 2 corresponds to the three layer cylinder 



3 and R = 1.0, 



< |d, < 0.3794 
0.3794 < |a;| < 0.5662 
0.5662 < |a;| < 0.6377 
0.6377 < \x\ < 1.0 

Let the data g(9) be collected for just one wave number fco = 3.0. Figures 1 and 2 
show the real and imaginary parts of the solutions for these two configurations. The 
solutions are practically indistinguishable, especially if noise is present. Letting Q% 
to be the original configuration for which the data g{9) is observed, the value of $ 
at the configuration Q 2 is just 0.00012. Thus, there is no way (by any method) to 
determine the original configuration Qi of the scatterer. Clearly, there are many 
more configurations that would produce practically identical observations. Even 
if it could be proven that, theoretically, there is a unique solution for this IPS, it 
would be useless in practice, because of this and other practically undistinguishablc 
configurations. 

On the other hand, the situation is quite different if we allow the scatterer to be 
probed with waves of multiple frequencies. 

Figures 3 and 4 show the real and imaginary parts for the same configurations Qi 
and Q 2 when the free space wave number ko is equal to 10.0. Then &(Q 2 ) = 1.4307. 
It is, of course, possible that there are configurations undistinguishablc at this 
frequency, but, combining the output for several frequencies we can hope to achieve 
a reasonable recovery of the original scatterer. We show in the subsequent sections, 
that it is, indeed, the case. While there are many theoretical questions concerning 



n(x) 



0.040 
8.282 
5.969 
1.0 
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the best, or a reasonable choice of frequencies, uniqueness for the IPM, stability 

estimates, etc., this work indicates the practicality of the multifrcquency approach. 

To illustrate this point further, let P be the set of three free space wave numbers 
(v) 

kg chosen to be 



P = {3.0, 6.5, 10.0}. (3.8) 

Figure 5 shows the profile of the functional $ as a function of the variable r , 0.1 < 
r < 0.6 in the configurations q r with 

< |x| < r 
r < \x\ < 0.6 
0.6 < |x| < 1.0 

The best fit to data functional exhibits a sharp minimum at r = .4, thus there is 
a hope to identify the sought configuration. 

4. LOCAL MINIMIZATION METHODS 

Using the best fit to data functional $ defined in (3.5), the IPM is reduced to 
a restrained minimization over the admissible set A a 4 m , defined in (3.6) and (3.7). 



n(x) 
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It is well known that a multidimensional minimization is an extremely difficult 
problem, unless the objective function is "well behaved". The most important 
quality of such a cooperative function is the presence of just a few local minima. 
Unfortunately, this is, decidedly, not the case in many applied problems, and, in 
particular, for the problem under the consideration. 

Figure 5 shows that our objective function $ has many local minima even along 
this arbitrarily chosen one dimensional cross-section of the admissible set. There 
are sharp peaks and large gradients. Consequently, the gradient based methods 
(see [7], [11], [14], [17], [19], [22]) would not be successful for a significant portion of 
this region. It is also appropriate to notice that the dependency of $ on its ar- 
guments is highly nonlinear. Thus, the gradient computations have to be done 
numerically, which makes them computationally expensive. More importantly, the 
gradient based minimization methods (as expected) perform poorly for these prob- 
lems. These complications are avoided by considering conjugate gradient type al- 
gorithms which do not require the knowledge of the derivatives at all. One such 
method is the Powell's method. 

For an TV-dimensional space the method can be described as follow (see [7]). 

Powell's Method 

1. Initialize the set of directions Ui to the basis vectors 



Ui = ej , i = 1,2, . ..,N 
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2. Save your starting position as Qq. 

3. For i = 1, . . . , N move to the minimum along the direction and call 
this point Qi. 

4. For i = 1, ... ,N, set Ui = 

5. Set un = Pn — Po- 

6. Move Pjv to the minimum along direction ujy and call this point Prj. 

It can be shown that an iteration of this procedure produces a set Uj of mutually 
conjugate directions, provided, as usual, that the objective function is quadratic. It 
also implies a quadratic convergence for nearly quadratic functions. The main dif- 
ficulty here is that the obtained set of conjugate directions tends to become "folded 
up", that is linearly dependent. However, as noted in [7], the set of directions Ui 
can be reset to the basis vectors ej after every N or N + 1 iterations of the basic 
procedure. 

As explained in the next Section, we leave the global exploration of the admissible 
set to global minimization methods. A local minimization is used to explore an 
immediate vicinity of the initial configuration Q £ K 2M . With this goal in mind, 
given a configuration Q e R 2M and a direction u in R 2M we seek a minimum of 
$ along this direction (by a bisection or a Golden Rule method) by restricting the 
probed points (at every minimization step) to the admissible set and by keeping 
them within a certain distance from the initial minimization point. This distance 
d m i n is determined a priori to be a percentage of the characteristic length of A a d m - 
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More precisely, the "turtle" one-dimensional minimization is done as follows. 
One- dimensional minimization 

1. Let the starting position be Qo- 

2. Move from Qo along the given direction u by the distance d m i n to obtain 
Qi G A a d m . 

3. Find the minimum of <J> on the interval [Qo, Q\]. 

If the minimum is attained inside the interval, then stop. 
If the minimum is attained at Qo, then reverse the direction. 

If the minimum is attained at Q\, then rename Qo — Qi, and repeat the procedure. 

We have used Brent's minimization Method [7] for the one-dimensional mini- 
mization in the step 3. This way the local minimum closest (the resolution is set 
up by d m in) to the starting configuration Q is determined. The choise of d min 
has to be balanced between the desire to explore the fine structure of the objective 
function and the computational costs. 

Now we can describe our Basic Local Minimization Method in K 2M . The above 
"turtle" onc-dimcnsional minimization procedure is used in all the minimization 
steps below. 

Basic Local Minimization Method 
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FIG. 5. Best fit profile for the configurations q r \ Multiple frequencies P = {3.0, 6.5, 10.0}. 

1. Initialize the set of directions ui to the basis vectors 

Ui = ei, i = 1,2,...,2M. 

2. Save your starting position as Qq. 

3. For i = 1, . . . , 2M move from Qq along the direction Ui to find the point of 
minimum Q\. 

4. Reindex the directions itj, so that (for the new indices) &(Q{) < ^{Qi) < 
,...,$(Q| M )<$(Qo). 

5. For i = 1, . . . , 2M move Qi-i to the minimum along the direction itj and call 
this point Qi. 

6. Set v = Q 2 m - Qq- 

7. Move Q2M to the minimum along direction v and call this point Qo- 

8. Repeat the above steps till a stopping criterion is satisfied. 

Note, that we use the temporary points of minima Q\ only to rearrange the 
initial directions itj in a different order. This method falls within the category of 
the Powell's minimization methods, and, as mentioned above, produces conjugate 
directions and a quadratic convergence for nearly quadratic functions. 
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Still another refinement of the above algorithm has turned out to be necessary 
to produce a successful minimization. Since the dimension 2M of the minimization 
space was chosen a priori to be larger than 27V, where TV is the (unknown) number 
of layers in the original scatterer, we expect that the sought point of minimum 
will be located in a lower dimensional subspace of the minimization space R 2M . 
This information available from the specific structure of our minimization problem 
appears to be nontrivial. Suffices to say, that all of our numerical experiments de- 
scribed in Section 6 have failed without the following (space dimension) "reduction" 
procedure. The main idea behind it is to conduct the local minimization searches 
in as low-dimensional subspaces as possible. It is specific to the inverse scattering 
problem for multilayer scatterer. 

If two adjacent layers have close refraction coefficients in the sense that the 
objective functional <& is not changed much when the two layers assigned the same 
refraction coefficient, then these two layers can be replaced with just one occupying 
their place. The minimization problem becomes constrained to a lower dimensional 
subspace of M 2M and the local minimization is done in this subspace. A similar 
procedure was used by us in [15] for the search of small subsurface objects. 

Reduction Procedure 

Let e r be a positive number. 

1. Save your starting configuration Q — (ri, r 2 , . . . , Tm, ni, «2, ■ • • , wm) and 
$(Qo). Let the M + 1-st layer be Dm+i — {cm < \x\ < R} and um+i — k^. 

2. For i = 2, . . . , M + 1 replace rn-\ in the layer Z?i_i by n,. Compute <E> at the 
new configuration Qf, and the difference cf = |$(Qo) — &(Qf)\- 

3. For i = 1, . . . , M replace n i+ i in the layer D i+ i by nj. Compute $ at the new 
configuration Qf, and the difference c" = |$(Qo) - ®(Qi)\- 

4. Find the smallest among the numbers cf and c". If this number is less than 
e r §(Qo), then adjust the refraction coefficient to rij in the "down" or "up" layer 
accordingly. Replace the two adjacent layers with one occupying their place, and 
renumber the layers. 

5. Repeat the above steps till no further reduction in the number of layers is 
occurring. 

Note, that an application of the Reduction Procedure may or may not result in 
the actual reduction of layers. 

Finally, the entire Local Minimization Method (LMM) consists of the following: 

Local Minimization Method (LMM) 

1. Let your starting configuration be Qo = (ri, T2, ■ ■ ■ , pmj n i, n 2, ■ ■ ■ , tim)- 

2. Apply the Reduction Procedure to Qo, and obtain a reduced configuration Qq 
containing M r layers. 

3. Apply the Basic Minimization Method in A a( im 

P| R 2Mr with the starting point 

Qq, and obtain a configuration Q\. 

4. Apply the Reduction Procedure to Q\, and obtain a final reduced configuration 

Q r v 
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5. GLOBAL MINIMIZATION METHODS 

Given an initial configuration Qo a local minimization method finds a local min- 
imum near Qo- On the other hand, global minimization methods explore the entire 
admissible set to find a global minimum of the objective function. While the lo- 
cal minimization is, usually, deterministic, the majority of the global methods are 
probabilistic in their nature. There is a great interest and activity in the devel- 
opment of efficient global minimization methods, see e.g. [4], [6]. Among them are 
the simulated annealing method ([20], [21]), various genetic algorithms [16], interval 
method, TRUST method ([2], [3]), etc. As we have already mentioned before, the 
best fit to data functional <& has many narrow local minima. In this situation it is 
exceedingly unlikely to get the minima points by chance alone. Thus our special 
interest is for the minimization methods, which combine a global search with a 
local minimization. In [15] we developed such a method (the Hybrid Stochastic- 
Deterministic Method), and applied it for the identification of small subsurface 
particles, provided a set of surface measurements. The HSD method could be clas- 
sified as a variation of a genetic algorithm with a local search with reduction. In 
this paper we consider the performance of two algorithms: Deep's Method, and 
Rinnooy Kan and Timmer's Multilevel Single-Linkage Method. Both combine a 
global and a local search to determine a global minimum. Recently these methods 
have been applied to a similar problem of the identification of particles from their 
light scattering characteristics in [29]. Unlike [29], our experience shows that Deep's 
method has failed consistently for the type of problems we are considering. See [10] 
and [29] for more details on Deep's Method. 

Multilevel Single-Linkage Method (MSLM) 

Rinnooy Kan and Timmcr in [25] and [26] give a detailed description of this 
algorithm. Zakovic et. al. in [29] describe in detail an experience of its application 
to an inverse light scattering problem. They also discuss different stopping criteria 
for the MSLM. Thus, we only give here a shortened and an informal description of 
this method and of its algorithm. 

In a pure Random Search method a batch H of L trial points is generated in 
A a dm using a uniformly distributed random variable. Then a local search is started 
from each of these L points. A local minimum with the smallest value of $ is 
declared to be the global one. 

A refinement of the Random Search is the Reduced Sample Random Search 
method. Here we use only a certain fixed fraction 7 < 1 of the original batch of L 
points to proceed with the local searches. This reduced sample H re d of jL points 
is chosen to contain the points with the smallest jL values of $ among the original 
batch. The local searches are started from the points in this reduced sample. 

Since the local searches dominate the computational costs, we would like to ini- 
tiate them only when it is truly necessary. Given a critical distance d we define a 
cluster to be a group of points located within the distance d of each other. Intu- 
itively, a local search started from the points within a cluster should result in the 
same local minimum, and, therefore, should be initiated only once in each cluster. 

Having tried all the points in the reduced sample we have an information on the 
number of local searches performed and the number of local minima found. This 
information and the critical distance d can be used to determine a statistical level of 
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confidence, that all the local minima have been found. The algorithm is terminated 
(a stopping criterion is satisfied) if an a priori level of confidence is reached. 

If, however, the stopping criterion is not satisfied, we perform another iteration 
of the MSLM by generating another batch of L trial points. Then it is combined 
with the previously generated batches to obtain an enlarged batch W of jL points 
(at iteration j), which leads to a reduced sample H 3 red of jjL points. The critical 
distance d is reduced to dj, (thus, the cluster's size is redefined), a local minimiza- 
tion is attempted once within each cluster, the information on the number of local 
minimizations performed and the local minima found is used to determine if the 
algorithm should be terminated, etc. 

The following is an adaptation of the MSLM method to the inverse scattering 
problem presented in Sections 2 and 3, with all the relevant notations. The LMM 
local minimization method introduced in the previous Section is used here to per- 
form local searches. 

MSLM 

(at iteration j). 

1. Generate another batch of L trial points (configurations) from a random uni- 
form distribution in A a( i m - Combine it with the previously generated batches to 
obtain an enlarged batch H 3 of jL points. 

2. Reduce H 3 ' to the reduced sample H J red of -yjL points, by selecting the points 
with the smallest 7j'L values of <& in H 3 . 

3. Calculate the critical distance dj by 



^ ^(v(^^ M \(, „ ]M <y\njL \ 1/M 

«j ■ = 7T ' I 1 I 1 + — \ [n hlgh - n low ) — I 

4. Order the sample points in H red so that < 3>(<3i+i) , i = 1,2,..., jjL. 
For each value of i, start the local minimization from Qi, unless there exists an 
index k < i, such that \\Qh — Qi\\ < dj. Ascertain if the result is a known local 
minimum. 

5. Let K be the number of local minimizations performed, and W be the number 
of different local minima found. Let 

w(K - l) 



K-W-2 
The algorithm is terminated if 



W tot < W + 0.5 . 



(5.1) 
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Here T is the gamma function, and a is a fixed constant. 

A related algorithm (the Mode Analysis) is based on a subdivision of the admis- 
sible set into smaller volumes associated with local minima. This algorithm is also 
discussed in [25] and [26]. From the numerical studies presented there, the authors 
deduce their preference for the MSLM. 

6. NUMERICAL RESULTS 

Introducing polar coordinates in (2.3) and separating the variables, equations for 
the total field u(x) become 

oo 

1— — 00 

for x 6 .Di, 



oo 

u m {x) = Y ( a m,iJi(k m \x\) + b m jYi(k m \x\))e 

l— — oo 

for x e Di, I = 2, . . . , N, and 



u{x) = e ik<x ' v> + AiH[ 1] {k \x\)e m 

l— — oo 

for x e D : tm < \x\ < R. Here J;,Y/ are the Bessel functions of the first and 
second kind, is the Hankel function of the first kind, and v is the direction 

vector of the incident wave. Since 

oo 
l — — oo 

for v = (1,0), the above equations and the conditions of continuity form a system 
of equations from which the field u{x) can be calculated on the circle S. This 
solves the direct problem (2.2)-(2.3). Other methods of solution for such problems 
are known as well, see e.g. [18], [28]. Solving the direct problem for the set P of 
three free wave numbers k^ (see (3.8)), we obtain the total fields u^(x). Their 
restrictions to S give the (simulated) data g( p \6). 

Our approach to the inverse problem IPM (see Section 2) is to recast it in the 
best fit to data form (3.5), and to minimize the objective functional $ over A a d m - 
We have tested Deep's global minimization method, the Multilevel Single-Linkage 
Method, and a Reduced Sample Random Search method. Each method was tried 
for three different original configurations Qo described below. The data g^ p \0) 
was computed at 120 angles 9i — 2ttI/120, I — 1,2,..., 120, and $ was evaluated 
according to (3.4) and (3.5). This data was used with three different noise levels 
S = 0.00,0.03 and 0.10. More precisely, for the uniformly distributed on [0,1] 
random variable z 
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gs(0)=g(6) + S\\g\\(2z-l)(l + i) 

for the noise level 6. 

Since our goal was to test the applicable algorithms, the values for the refrac- 
tion coefficients, the size, the wave numbers, etc., were chosen arbitrarily at this 
time, that is without a regard for their possible physical relevancy. The original 
configurations are: 

Configuration . This is a one layer cylinder = (0.72, 4.2025) with N = 1 
and R — 1.0, see (3.2), that is the refraction coefficient is defined by 

[4.2025 < \x\ < 0.72 

n \ x ) = \ 

[1.0 0.72 < \x\ < 1.0 

Configuration Q c ( , 2) . This is a two layer cylinder Qq 2) = (0.4,0.6,0.49,9.0) with 
N = 2 and R = 1.0, that is 

< \x\ < 0.4 
0.4 < |a;| < 0.6 
0.6 < \x\ < 1.0 

Configuration Q< 3) . Three layer cylinder = (0.3, 0.7, 0.8, 4.0, 25.0, 9.0) with 
N = 3 and R = 1.0, that is 

< \x\ < 0.3 
0.3 < |a;| < 0.7 
0.7 < \x\ < 0.8 
0.8 < \x\ < 1.0 

To identify these configurations we applied the global minimization methods of 
Section 5. In each one we let M — 4 , R = 1.0. A priori bounds for the refraction 
coefficients were chosen to be ni ow = .04 and rihigh — 30.25. Minor modifications 
to the description of the methods in Section 5 were introduced for the purpose of 
computational simplification. In particular, the minimization was done in ^Jn m 
rather than in n m as stated there. This results in a rescaling of the admissible 
set. In each case, after a global minimum Qjmin was determined, the error of the 
identification 



n{x) 




n(x) 



4.0 
25.0 
9.0 
1.0 



_ J D \n mm (x)-n(x)\ 

J D n(x) ' (bA) 

where n(x) is the refraction coefficient of the original configuration Q , was com- 
puted to determine if the identification was successful. We distinguished between 
the two levels of a successful identification: e err < 0.01 and e err < 0.1. 
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Identification by Deep's Method([10j, [29]) 

Each test of the method consisted in 100 independent runs. Since M = 4 the min- 
imization was done in M 8 . As we have already mentioned, the method failed every 
time. It seems, that the local minimization phase of Deep's method (minimization 
over randomly selected parabolas) is not extensive enough to identify narrow local 
minima present in this problem. Also, the method docs not use the Reduction 
Procedure (see Section 4), which, we think, is another reason for its failure. As in 
[29] we have also observed the cycling of the algorithm. 

Identification by Reduced Sample Random Search Method 

This method is presented in Section 5 in the subsection on the Multilevel Single- 
Linkage Method. The Local Minimization Method with the Reduction Procedure 
(as described in Section 4) was used in the local minimization phase. Chosen pa- 
rameters L — 15000 and 7 = 0.01 the performance of this method is the same as the 
Multilevel Single-Linkage Method. In fact, L = 15000 is exactly the sample size in 
MSLM at it termination in our experiments. Since MSLM has the great advantage 
of a self-contained statistical stopping criteria (and from which the number 15000 
was determined in the first place), it is, clearly, a preferred method. 

Identification by Multilevel Single-Linkage Method 

We have attempted to identify all 3 original configurations Qq~\ and . 
Each one with no noise in the data (6 = 0.00) as well as with noise levels 8 = 0.03 
and 5 = 0.10. Each of the 9 tests consisted of 10 independent runs. It took about 
60 to 80 minutes on average to complete one run on a 333 MHZ PC. We used 
M = 4 , R = 1.0 , 7 = 0.01 and the sample size L = 200. The parameter a was 
chosen to be equal to 1.0. Value a = 4.0 was used in [26], and a = 1.9 in [12]. As 
in Deep's Method above, a priori bounds for the refraction coefficients were chosen 
to be ni ow — .04 and rihigh — 30.25. The value e r = 0.1 was used in the Reduction 
Procedure (see Section 4) during the local minimization phase. 

As in other works on the clustering algorithm, we have found the stopping rule 
(5.1) to be unsatisfactory. In our experience the difference Wtot — W, while slightly 
decreasing with the number of performed minimizations, has quickly stabilized 
around the value of 5. Thus, the stopping criterion (5.1) 

Wtot < W + 0.5 . 

could not be attained. This issue has been discussed in [12] and [5], where 
a different stopping rule was suggested for functions with large number of local 
minima. Since the Bayesian stopping rule reflects the level of confidence in finding 
all the local minima, a relaxation of (5.1) would mean a smaller level of confidence, 
which still may be acceptable to assure that the global minimum is found among 
already performed local minimizations. We have chosen to replace (5.1) with 

Wtot <W + 0.5 or Wtot < (l + e tot )W, (6.2) 

where e to t — 0.03. 
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TABLE 1 

Identification by MSLM. Original configuration Qq 1 '. 

Success rate 

Noise Runs 0.01 0.1 Smallest $ 

J=0.00 10 8 10 0.0000 

5=0.03 10 9 10 0.0016 

(5=0.10 10 10 10 0.0178 



TABLE 2 

Identification by MSLM. Original configuration . 



Noise 


Runs 


Success rate 1 

0.01 


0.1 


Smallest $ 


(5=0.00 


10 


10 


10 


0.0000 


(5=0.03 


10 


1 


10 


0.0034 


(5=0.10 


10 


2 


9 


0.0362 



As before 



K-W-2' 



where K is the number of local minimizations performed, and W is the number of 
different local minima found. Thus, the MSLM algorithm is terminated if (6.2) is 
satisfied. In our numerical experiments we have got the following average values 
K = 5000 , W = 150, and W tot = 155. 

(2) 

An example of a successful (e err < 0.1) identification for Q y ' and 5 = 0.10 is 
shown on Figure 6. The identified configuration is a two layer cylinder 



Q ld = (0.3966, 0.5943, 0.4684, 9.203) , 
with <&{Q ld ) = 0.0367, e err = 0.0480. That is 



n(x) = < 



0.4684 < \x\ < 0.3966 
9.203 0.3966 < \x\ < 0.5943 
1.0 0.5943 < Id < 1.0 



An example of a successful (e err < 0.1) identification for and 5 = 0.03 is a 
three layer cylinder 



1 Identification is successful if e err < 0.01, or e err < 0.1 correspondingly, see (6.1). 
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r 

(2) 

FIG. 6. Refraction coefficients n(x) for the original Q () and the identified Qid (solid line) 
configurations. Data noise level <S = 0.10, *(Q id ) = 0.0367, e err = 0.0480. 



TABLE 3 

Identification by MSLM. Original configuration Qq . 







Success rate 1 




Noise 


Runs 


0.01 0.2 


Smallest $ 


5=0.00 


10 


2 7 


0.0000 


5=0.03 


10 


5 


0.0071 


5=0.10 


10 


5 


0.0541 




Q ld = (0.3030, 0.7067, 0.8079, 4.071, 24.528, 8.857) , 




with 


i) = 0.00708, e err 
n(x) = < 


= 0.04125. That is 

'4.071 < \x\ < 0.3030 
24.528 0.3030 < \x\ < 0.7067 
8.857 0.7067 < |a;| < 0.8079 
1.0 0.8079 < \x\ < 1.0 
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7. CONCLUSIONS 

The inverse scattering problem IPM is the identification of a multilayered scat- 
terer by a set of observations on its boundary. Such problems have applications in 
science and engineering. In the case of the weak scattering approximation, many 
such problems can be solved by a linearized inversion. However, if the scattering is 
not weak, other methods of solution need to be developed. We have illustrated in 
Section 3, that an inversion based on just one frequency of the incident waves can- 
not be successful, since there are distinct configurations, producing practically the 
same observations. Introducing multiple frequencies, however, makes the inverse 
problem more amenable to a solution. 

In this paper the inverse problem is transformed into the best fit to data min- 
imization problem. This minimization is difficult, since the objective function is 
rugged and has many narrow local minima. A promising way to treat such a min- 
imization is by a combination of global (probabilistic) and local (deterministic) 
minimization methods. In this paper we examined various local and global meth- 
ods. Concerning the local minimization methods it was shown, that the Local 
Minimization Method (LMM) of Section 4 was successful, even where other consid- 
ered methods failed. This method is a variation of a conjugate directions method 
with no use of partial derivatives. It has a quadratic convergence near quadratically 
shaped minima. However, even this method needs to be enhanced by a Reduction 
Procedure (Section 4) . This procedure helps the minimization to take an advantage 
of the a priori available information, that the sought minima are likely to be found 
in certain lower dimensional subspaces of the entire minimization space. 

For the global minimization part we considered Deep's method and the Multilevel 
Single-Linkage Method. While Deep's method failed, the MSLM was successful in 
many instances. It also has an important advantage of having termination criteria 
establishing a level of confidence, that the found minima contain the sought global 
minimum. Among the deficiencies of the MSLM are its slow execution, and incon- 
sistency and failturc to identify some configurations. There is still a problem in 
choosing an appropriate stopping rule. Thus, the MSLM provides a benchmark, 
against which the performance other methods can be judged and measured. 
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